% this version (update from v2.1) is set up for use with the looping
% function of scanimage, so for a single experiment with 2 repeats of 3
% stimuli, for example, you have 6 separate movies from scanimage + 1
% metadata file + 1 eye tracking movie + a 6x1 cell array for the ball
% tracking data


warning('off','all')
dostabilize = false;             %adds a image stabilization step if necessary
deinterleave = 2;               %if deinterleave == 1 then deinterleave the movies and only use the first channel, if deinterleave == 2 then use second channel
objtype = '16x_v2';        %'25x' or '16x' or '16x_v2' or '25x_v2'

if 1      %can skip the selection code if you've already done it
% dependencies:
%
%% few ways to select appropriate movies:

movDir = uigetdir(cd,'select directory containing the movies');

%deinterleave = 0 means only one channel collected;
%deinterleave = 1 means two channels collected but use channel 1 for
%tracking
%deinterleave = 2 means two channels collected but use channel 2 for
%tracking

   
cd(movDir)

fn = uigetfile('*.tif','select appropriate movies','multiselect','on');     %select only one movie in the sequence
if ischar(fn)
    temp = fn;
    fn = cell(1);
    fn{1} = temp;
end
filenames = struct;
for n = 1:length(fn)
    filenames(n).name = fn{n};
end

%% 
if 1    %preselect vessel centers
    for n = 1:length(filenames)
        movname = filenames(n).name;
        savename = strrep(movname,'.tif','_manualEndpointSelection.mat');
        if isfile(savename)     %this prevents re-doing existing selections
            continue;
        end
        try
            meta = getSImetadata(movname);
            curzoom = meta.zoom;
        catch
            disp('issue with metadata')
            continue;
        end
        if length(meta.savedchannels) > 1 && deinterleave == 0
            curinterleave = 1;
        elseif length(meta.savedchannels) == 1
            curinterleave = 0;
        else
            curinterleave = deinterleave;
        end
        switch curinterleave
            case 0
                try
                im = mean(loadTiffStack2_SI(movname,31:40),3);
                catch
                    'a'
                end
            case 1
                im = loadTiffStack2_SI(movname,31:40);
                im = mean(im(:,:,:,1),3);
            case 2
                temp = loadTiffStack2_SI(movname,31:40);
                im = mean(temp(:,:,:,2),3);
                refim2 = mean(temp(:,:,:,1),3);
                
                figure; imagesc(refim2);title('hydrazide reference');axis equal
        end
        
        disp(['using this objective settings: ',objtype])
        micronsPerPixel = pixel2micron_2p(curzoom,objtype);
        [ep1,ep2,endpointIDs] = linescansFromCentersMalabSelect(im,micronsPerPixel,movname);        %see comment in this function for changing the length of the linescan
        
        savename = strrep(movname,'.tif','_manualEndpointSelection.mat');
        save(savename,'ep1','ep2','endpointIDs','movname')
    end
end
    

%% port the vessel selections from the first technical repeat to the other repeats
% this is designed for the whisker stim paradigm but in principle could be
% useful for visual as well...
copyVesselSelectionsAcrossTrials;
else
    movDir = cd;
end
%% iterate through movies and track
posthocAveragingGlobal = 1;        % if you want to do additional time averaging (e.g. if you accidentally had averaging set to 1 instead of 3)
aquisitionAveraging = [];    % leave this as empty if you want to use the metadata frame averaging or set it to a value if you know the frame averaging beforehand that the software is doing
filenames = dir('*_manualEndpointSelection.mat');
for n = 1:length(filenames)
    temp = strrep(filenames(n).name,'_manualEndpointSelection.mat','_tracerTracking_preselected2.mat');
    if isfile(temp)
        load(temp);
        if isfield(results,'kymos') && ~isempty(results.kymos)      %don't re-do tracking if it already exists
            clear results
            continue;
        end
    end
    movname = strrep(filenames(n).name,'_manualEndpointSelection.mat','.tif');
    try
        [mov,meta] = loadTiffStack2_SI(movname,[],aquisitionAveraging);           %% note this 
    catch
        mov = loadTiffStack2(movname);
        meta = getSImetadata(movname);
    end
    
    if meta.frameaveraging == 1 && round(meta.framerate,0) == 30
        posthocAveraging = 3;
    else
        posthocAveraging = posthocAveragingGlobal;
    end
        
    
    disp(['using settings for this objective: ',objtype])
    pix2mic = pixel2micron_2p(meta.zoom,objtype);
    if length(meta.savedchannels) > 1 && deinterleave == 0
        curinterleave = 1;
    elseif length(meta.savedchannels) == 1
        curinterleave = 0;
    else
        curinterleave = deinterleave;
    end
    
    try
        load(strrep(movname,'.tif','_manualEndpointSelection.mat'));
        ind = find(~isnan(ep1(:,1)));
        ep1 = ep1(ind,:); ep2 = ep2(ind,:);endpointIDs = endpointIDs(ind);
        endpointStruct.ep1 = ep1;
        endpointStruct.ep2 = ep2;
        endpointStruct.endpointIDs = endpointIDs;
        if isempty(ep1)
            disp('no vessels selected for this one')
            results = [];
            continue;
        end
    catch
        disp('no vessels selected for this one')
        results = [];
        continue;
    end
    
    switch curinterleave
        case 1
            if meta.si_version == 2022
                mov = mov(:,:,:,1);
            else
                mov = mov(:,:,1:2:end);     %uses red channel
            end
        case 2
            if meta.si_version == 2022
                mov = mov(:,:,:,2);
            else
                mov = mov(:,:,2:2:end);     %uses green channel
            end
    end
    if dostabilize
        refframe = 1;
        mov = driftCorrectWidefield(mov,refframe);
        disp('drift correction done')
    end
    
    if posthocAveraging > 1
        mov = movmean(mov,posthocAveraging,3);
        mov = mov(:,:,(ceil(posthocAveraging/2):posthocAveraging:end));
        fps = meta.framerate/posthocAveraging;
    end
    
    try

        results = tracerTrackRecipe_v3_3(mov,endpointStruct,0,pix2mic);
        
    catch
        disp('no preloaded results')
        results = [];
        continue;
    end
    results.movname = movname;
    results.directory = movDir;
    results.objtype = objtype;
    %% go through each vessel and get its position in the whole coordinate system:
    curid = results.vesselID;
    uid = unique(curid);
    vesselPositions = zeros(length(uid),3);
    movmeta = getSImetadata(results.movname);
    vesselPositions(:,3) = movmeta.position(3);
    disp(['using settings for this objective: ',objtype])
    pix2mic = pixel2micron_2p(meta.zoom,objtype);
    for k = 1:length(uid)
        ind = find(results.vesselID == uid(k));
        ep = results.endpoints(ind,:,:);
        ep = mean(ep,3);                %calculate midline position of the vessel
        ep = ep*pix2mic;                %convert to microns
        offset = [movmeta.Height/2,movmeta.Width/2]*pix2mic;
        ep = ep-offset;
        
        temp = ceil(size(ep,1)/2);        %take middle of the tracked segment
        vesselPositions(k,1:2) = ep(temp,:) + movmeta.position(1:2);
    end
    results.vesselPositions = vesselPositions;
    
    %% save the results structure
    savename = strrep(movname,'.tif','_tracerTracking_preselected2.mat');
    save(savename,'results');
    clear results
end

warning('on','all')
